R Packages

library(sp)
library(gstat)
library(sf)
library(rgdal)
library(ggplot2)
library(plotly)
library(Matrix)

Data.

Air Poultion data CDMX 2020/01/16 17

datos <- read.csv("Air_polution_cdmx_2020_01_16_17h.csv")
datos <- datos[c("Estacion",
               "X",
               "Y",
               "NO2",
               "O3",
               "NOX")]

pander::pander((datos))
Estacion X Y NO2 O3 NOX
AJU 482901 2117907 NA 50 NA
AJM 478188 2130946 5 51 8
ATI 473346 2164689 20 70 20
CAM 482180 2152665 23 83 24
CCA 481502 2136931 6 46 8
CUA 469366 2141275 9 56 10
CUT 479189 2180751 13 75 14
FAC 474444 2154232 31 71 35
HGM 484020 2146380 25 81 31
IZT 487647 2143367 21 61 23
LLA 495842 2164872 14 81 15
LPR 487650 2160000 NA 64 NA
MER 487445 2147815 23 79 26
MGH 478716 2145543 NA 66 NA
MON 510196 2151776 9 80 10
MPA 498809 2123036 1 50 NA
NEZ 497038 2144394 8 62 9
PED 478557 2136817 6 56 6
SAG 496819 2159801 16 80 16
SFE 472393 2140390 9 66 10
TAH 498890 2128098 3 47 3
TLA 478535 2159383 26 64 30
TLI 481421 2167509 18 86 20
UAX 489113 2134517 NA 54 NA
UIZ 492241 2140751 7 45 7
VIF 489875 2173664 19 80 22
XAL 491355 2159031 27 80 27

Corregionlization

Semipositive definite matrix spherical model

mat1 <- cbind(c(30, 30, 30),
              c(30, 50, 30),
              c(30, 30, 35))
#matriz definida positiva "cercana"
mat1 <- data.frame(as.matrix(nearPD(mat1)$mat))
names(mat1) <- c("NO2", "O3", "NOX")
row.names(mat1) <- c("NO2", "O3", "NOX")
pander::pander(mat1)
  NO2 O3 NOX
NO2 30 30 30
O3 30 50 30
NOX 30 30 35

Semipositive definite matrix spherical hol.

mat2 <- cbind(c(13.02, 24.5, 18.739),
              c(24.58, 46.4, 35.36),
              c(18.73, 35.36, 26.95))
mat2 <- data.frame(as.matrix(nearPD(mat2)$mat))
names(mat2) <- c("NO2", "O3", "NOX")
row.names(mat2) <- c("NO2", "O3", "NOX")
pander::pander(mat2)
  NO2 O3 NOX
NO2 13.02 24.54 18.73
O3 24.54 46.4 35.36
NOX 18.73 35.36 26.96

G_stat object

Semivariogram

vgmno2 <- vgm(psill = mat1[1, 1],
            model = "Sph",
            range = 6096,
            add.to = vgm(psill = mat2[1, 1],
                         model = "Hol",
                         range = 2294))

vgmo3 <- vgm(psill = mat1[2, 2],
            model = "Sph",
            range = 6096,
            add.to = vgm(psill = mat2[2, 2],
                         model = "Hol",
                         range = 2294))

vgmnox <- vgm(psill = mat1[3, 3],
            model = "Sph",
            range = 6096,
            add.to = vgm(psill = mat2[3, 3],
                         model = "Hol",
                         range = 2294))

Cross semivarogram

vgmno2_o3 <- vgm(psill = mat1[1, 2], model = "Sph",
            range = 6096,
            add.to = vgm(psill = mat2[1, 2],
                         model = "Hol",
                         range = 2294))

vgmno2_nox <- vgm(psill = mat1[1, 3],
            model = "Sph",
            range = 6096,
            add.to = vgm(psill = mat2[1, 3],
                         model = "Hol",
                         range = 2294))

vgmno3_nox <- vgm(psill = mat1[2, 3],
                  model = "Sph",
                  range = 6096,
                  add.to = vgm(psill = mat2[2, 3],
                               model = "Hol",
                               range = 2294))

Gstat

remove_na <- function(frame, vari_) {

    # Remove na from sp object

    datos1 <- frame

    bool <- !is.na(datos1@data[vari_])
    datos1@data <- datos1@data[bool, ]
    datos1@coords <- datos1@coords[bool, ]

    return(datos1)

}

coordinates(datos) <- ~ X + Y

g_st <- gstat(NULL,
              id = "NO2",
              formula = NO2 ~ X + Y,
              model = vgmno2,
              data = remove_na(datos, "NO2"))

g_st <- gstat(g_st,
              id = "O3",
              formula = O3 ~ Y,
              model = vgmo3,
              data = remove_na(datos, "O3"))

g_st <- gstat(g_st,
              id = "NOX",
              formula = NOX ~ Y,
              model = vgmnox,
              data = remove_na(datos, "NOX"))
#cross


g_st <- gstat(g_st,
              id = c("NO2", "O3"),
              model = vgmno2_o3)

g_st <- gstat(g_st,
              id = c("NO2", "NOX"),
              model = vgmno2_nox)

g_st <- gstat(g_st,
              id = c("O3", "NOX"),
              model = vgmno3_nox)


pander::pander(do.call(rbind, g_st$model)[, 1:3])
  model psill range
NO2.1 Hol 13.02 2294
NO2.2 Sph 30 6096
O3.1 Hol 46.4 2294
O3.2 Sph 50 6096
NOX.1 Hol 26.96 2294
NOX.2 Sph 35 6096
NO2.O3.1 Hol 24.54 2294
NO2.O3.2 Sph 30 6096
NO2.NOX.1 Hol 18.73 2294
NO2.NOX.2 Sph 30 6096
O3.NOX.1 Hol 35.36 2294
O3.NOX.2 Sph 30 6096

Plot model model and estimation

plot(variogram(g_st),
     model = g_st$model,
     pl = T)

Prediction plot.

prediction_plot <- function(g_object, variable, map_path) {

    map <- readOGR(map_path)
    new <- sp::spsample(map, n = 100000, type = "regular")
    coordinates(new) ~ x1 + x2
    colnames(new@coords) <- c("X", "Y")

    predic <- predict(g_object, newdata = new)

    prediction <- data.frame(predic)

    pred <- paste(variable, ".pred", sep = "")

    plot <- ggplot(prediction, aes_string("X", "Y", fill = pred)) +
            geom_tile() +
            scale_fill_viridis_c() +
            theme_void()

    return(plot)

}


variance_plot <- function(g_object, variable, map_path) {

    map <- readOGR(map_path)
    new <- sp::spsample(map, n = 10000, type = "regular")
    coordinates(new) ~ x1 + x2
    colnames(new@coords) <- c("X", "Y")

    predic <- predict(g_object, newdata = new)

    prediction <- data.frame(predic)

    var <- paste(variable, ".var", sep = "")

    plot <- ggplot(prediction, aes_string("X", "Y", fill = var)) +
            geom_tile() +
            scale_fill_viridis_c(option = "inferno",
                                 direction = -1) +
            theme_void()

    return(plot)

}

cv_plot <- function(g_object, variable, map_path) {

    map <- readOGR(map_path)
    new <- sp::spsample(map, n = 10000, type = "regular")
    coordinates(new) ~ x1 + x2
    colnames(new@coords) <- c("X", "Y")

    predic <- predict(g_object, newdata = new)

    prediction <- data.frame(predic)
    pred <- paste(variable, ".pred", sep = "")
    var <- paste(variable, ".var", sep = "")
    aux <- abs(sqrt(prediction[var]) / abs(prediction[pred]))
    aux[aux > 1] <- 1
    prediction["cv"] <- aux

    plot <- ggplot(prediction, aes_string("X", "Y", fill = "cv")) +
            geom_tile() +
            scale_fill_viridis_c(option = "magma",
                                 direction = -1) +
            theme_void()

    return(plot)
}


pl1 <- prediction_plot(g_st, "O3",
                       "SP/mpiosutm.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/juan/Documents/Inves/Webinar/Day1/2_COK_G_stat/SP/mpiosutm.shp", layer: "mpiosutm"
## with 54 features
## It has 7 fields
## Linear Model of Coregionalization found. Good.
## [using universal cokriging]
pl2 <- variance_plot(g_st, "O3",
                     "SP/mpiosutm.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/juan/Documents/Inves/Webinar/Day1/2_COK_G_stat/SP/mpiosutm.shp", layer: "mpiosutm"
## with 54 features
## It has 7 fields
## Linear Model of Coregionalization found. Good.
## [using universal cokriging]
pl3 <- cv_plot(g_st, "O3",
               "SP/mpiosutm.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/juan/Documents/Inves/Webinar/Day1/2_COK_G_stat/SP/mpiosutm.shp", layer: "mpiosutm"
## with 54 features
## It has 7 fields
## Linear Model of Coregionalization found. Good.
## [using universal cokriging]
ggplotly(pl1)
ggplotly(pl2)
ggplotly(pl3)